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ABSTRACT 


Equivalent  impedance  of  a  rough  surface  is  determined  by  solving  the  problem  of 
radiowave  propagation  over  rough  random  surface  by  using  the  parabolic  equation  method. 
We  consider  horizontal  polarization  and  assume  the  rough  surface,  defined  by  discrete 
points,  to  be  perfectly  conducting.  This  assumption  has  minimum  effects  for  frequencies  in 
the  VHF  band  and  above.  The  essential  point  of  this  thesis  is  to  solve  numerically  a 
Volterra  integral  equation  of  the  second  kind  for  the  surface  current.  Because  the 
computation  becomes  very  intense  for  long  ranges,  we  present  a  multiple  sections  solution, 
where  the  solution  procedure  is  successively  repeated  for  smaller  ranges.  The  numerical 
results  are  compared  with  results  available  in  the  literature.  Finally,  we  apply  the  technique 
to  determine  the  equivalent  impedance  of  a  sinusoidal  random  rough  surface  by  comparing 
the  resulting  normalized  field  to  the  field  of  a  2-ray  model  over  a  constant  impedance  plane. 
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I.  INTRODUCTION 


Low-grazing  angle  radiowave  propagation  over  random  rough  surface  is  an 
extremely  important  topic  in  many  areas  of  interest  to  the  Navy,  such  as  sea-to-sea 
communications  or  detection  of  low-flying  targets  over  the  sea  surface.  Depending  on 
the  rough  surface  characteristics,  the  direct  field  and  the  field  reflected  and  diffracted  by 
the  rough  surface  will  interact  in  different  ways,  changing  the  field  strength  at  a  certain 
point  of  interest.  Numerous  analytical  methods  are  available  for  predicting 
electromagnetic  wave  propagation  over  such  surfaces,  such  as  physical  optics,  geometric 
optics,  normal  mode  analysis  and  combination  of  the  above.  However,  a  complex  surface 
profile  complicates  the  application  of  these  methods.  The  Parabolic  Equation  (PE) 
method,  where  a  Helmholtz  equation  is  approximated  by  a  parabolic  equation,  has  been 
extensively  used  for  the  determination  of  forward  propagation  because  of  its  numerical 
efficiency.  When  making  this  approximation  we  are  assuming  that  waves  travel 
predominantly  in  the  forward  direction  and  that  the  back-scattered  waves  are  neglected. 
This  situation  is  met  in  many  practical  propagation  problems.  Naturally  there  are 
situations  where  it  is  not  possible  to  apply  such  an  approximation.  For  example,  when 
there  is  a  large  obstacle  behind  a  receiver,  the  back-scattered  field  cannot  be  neglected. 
The  use  of  the  parabolic  partial  differential  equation  method  allows  us  to  calculate  the 
field  at  one  location  in  terms  of  the  field  at  a  previous  location,  utilizing  a  marching 
scheme,  resulting  in  a  more  efficient  computational  method. 

In  this  thesis  we  numerically  solve  a  Volterra  integral  equation  of  the  second  kind 
to  predict  fields  over  a  rough  surface.  We  consider  horizontal  polarization  and  a 
perfectly  conducting  rough  surface  to  simplify  the  formulation.  The  later  assumption  will 
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have  minimum  effects  for  frequencies  in  the  VHF  band  and  above.  Then  we  apply  the 
technique  developed  to  determine  the  equivalent  impedance  of  sinusoidal  rough  surfaces. 
This  is  done  by  comparing  the  resulting  normalized  field  to  the  field  of  a  2-ray  model 
over  a  constant  impedance  plane.  A  random  surface  with  a  fixed  spectrum  is  simulated 
by  translating  the  sinusoidal  surface  of  given  peak-to-peak  height  and  period  in  a  random 
fashion.  This  is  done  by  including  a  phase  for  the  sinusoidal  surface  and  varying  it 
uniformly  over  [0,2/r] .  Finally  we  show  how  to  generate  the  rough  sea  surface  according 
to  the  Pierson-Moskowitz  spectrum. 

Chapter  II  describes  the  theoretical  formulation,  the  numerical  solution  used  to 
solve  the  integral  equation  and  the  methodology  for  the  equivalent  impedance  calculation 
when  using  parabolic  equation.  Chapter  IE  presents  numerical  results,  comparing  them 
with  practical  cases  of  radiowave  propagation  over  terrain,  and  calculates  the  equivalent 
impedance  for  a  sinusoidal  surface.  Chapter  IV  describes  how  to  generate  the  sea  surface 
according  to  the  Pierson-Moskowitz  spectrum.  Recommendations  and  conclusions  are 
presented  in  Chapter  V. 
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II.  THEORETICAL  FORMULATION 


In  this  chapter  we  first  present  a  Volterra  type  integral  equation  of  the  second  kind 

for  the  determination  of  the  surface  current  on  the  rough  surface  based  on  the  parabolic 

equation  approximation.  Because  a  straightforward  solution  of  a  Volterra  integral 

equation  for  long  ranges  becomes  computationally  burdensome,  we  present  an  alternate 

technique  whereby  (i)  the  range  is  divided  into  several  sections,  (ii)  the  field  is  calculated 

at  vertical  lines  at  the  junction  of  two  sections,  and  (iii)  the  field  on  a  vertical  line  is  used 

as  a  source  for  the  determination  of  current  on  the  next  section.  The  rough  surface  is 

defined  by  discrete  points  and  a  linear  interpolation  between  consecutive  points  is 

assumed.  Figure  1  shows  a  possible  irregular  surface.  The  terrain  is  assumed  to  be  a 

perfect  electric  conductor  (conductivity  —>  infinite).  A  2-dimensional  case  is  considered 

where  the  source,  geometry  and  fields  are  assumed  to  be  invariant  in  the  y  direction. 

Propagation  is  assumed  to  take  place  in  the  xz  plane  and  is  mathematically  described  by 

the  standard  parabolic  equation  given  by 

du  _  i  d2u  .. 

~te~2k0  d?  {  ' 

where  V  is  the  reduced  field  variable  and  is  related  to  the  electrical  field  ’[/’  through 

U(x,z)  =  e,k°xu(x,z), 


is  the  wave  number  and  A  is  the  free-space  wavelength.  An  e  "“time 


dependence  is  assumed  throughout. 
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The  following  assumptions  are  made:  (i)  the  waves  are  assumed  to  predominantly 
travel  in  the  positive  x-direction,  and  (ii)  back  scattering  of  waves  is  ignored.  Given  a 
source  field,  f(z),  at  x  =  0  one  wants  to  calculate  the  field  at  x,  where  x  is  the 
horizontal  distance  from  the  source  and  z  is  the  height  with  respect  to  some  reference 
level.  An  integral  representation,  similar  to  Kirchoff’s  integral  form,  can  be  derived 
based  on  integrating  the  parabolic  equation  in  the  domain  formed  by  two  vertical  lilies, 
the  rough  surface,  and  an  infinite  boundary  as  shown  in  Figure  1. 


Figure  1 :  Perfectly  Conducting  Rough  Surface 
The  field  at  the  point  (x,  z)  is  given  by  [Ref.  1]: 

u(x,  z)  =  /  07 )G0  [*,  z;0,  T]]dT]  -  Jq  [x,z;^,JJ  =  g  (£)]d4 

( lef 

=  u1(x,z)-u2(x,z)  (2-2) 


4 


where  //(•)  is  the  Heaviside  step  function. 

The  first  integral  of  Equation  2.2  is  performed  in  the  vertical  direction  using  the 
source  field  f(rj) .  Its  lower  limit  is  g0  and  its  upper  limit  is  infinite.  However,  when  a 
numerical  computation  is  required,  it  will  be  truncated  at  some  finite  height.  This 
integral  corresponds  to  the  field  produced  by  the  source  f(tj)  in  free  space. 

The  second  integral  of  Equation  2.2  is  performed  in  the  horizontal  direction  along 
the  rough  surface  T]  —  g(£).  This  integral  corresponds  to  the  field  reflected  and 

du 

diffracted  by  the  rough  surface.  The  term  -—*-(£)  corresponds  to  the  vertical  derivative 

OTJ 

of  the  field  on  the  surface  and  is  related  to  the  normal  derivative  on  the  rough  surface  by 

1  3m  ^ 

where  g(x)  is  the  surface  profile,  and  the  prime  denotes  differentiation  with  respect  to 
the  argument.  Whenever  we  refer  to  the  term  current  we  will  assume  it  to  mean  the 
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du 

vertical  derivative  of  surface  field.  An  integral  equation  for  — —  may  be  formulated  by 

oij 

differentiating  Equation  2.2  with  respect  to  z  and  tending  z  to  the  surface.  Taking  the 
proper  limiting  values  one  gets  [Ref.  1]: 

dU°  (x)  =  -2  j  f(7])~G0[x,  g(x);0, T]]dTj+ 


Bz 


So 


+  f g(X) - &)-G,lx,g(xy,l  y  tf) Vi  (2.6) 

J0  Bri  x-£ 

where  denotes  a  principal  value  integral. 

Equation  2.6  is  a  Volterra  integral  equation  of  the  second  kind  where  the  current 
at  a  given  point  will  depend  only  on  the  current  at  the  previous  points.  Thus  we  notice 
the  causal  property  of  the  Volterra  integral  equation.  In  order  to  solve  it  we  need  to  use  a 
marching  procedure,  where  we  calculate  the  current  at  a  point  based  on  the  current  at  all 
other  previous  points.  However,  this  current  calculation  becomes  computationally  very 
intense  for  large  distances.  In  order  to  alleviate  this  problem,  we  present  the  multiple 
sections  solution  described  on  the  next  section.  Once  we  have  the  current  on  the  surface, 
given  by  Equation  2.6,  we  can  find  the  field  at  any  point  (x,  z)  using  Equation  2.2. 

A.  MULTIPLE  SECTIONS  SOLUTION 

We  can  reduce  the  computation  time  by  dividing  the  surface  into  a  number  of 
sections  and  using  Equation  2.6  over  the  sub-intervals.  Thus,  instead  of  calculating  the 
current  on  the  entire  surface  between  a  transmitter  and  a  location  of  interest,  we  can 
perform  current  calculations  successively  on  shorter  sections  by  making  use  of  field  data 
available  on  vertical  lines. 
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Given  the  field  on  a  vertical  line,  the  following  procedure  is  performed  in  order  to 
calculate  the  field  on  the  next  vertical  line: 

Step  1:  Calculate  the  incident  current  induced  on  the  rough  surface  by  the  initial 
field.  The  incident  current  on  the  surface  is  given  by  the  first  term  on  the  right  hand  side 
in  Equation  2.6. 

Step  2:  Calculate  the  total  current  on  the  rough  surface.  The  total  current  on  the 
surface  is  given  by  both  terms  on  the  right  hand  side  of  Equation  2.6. 

Step  3:  Calculate  th e  field  at  the  next  vertical  line,  utilizing  Equation  2.2. 

We  notice  a  trade  off  when  applying  this  split-section  method,  illustrated  on 
Figure  2.  On  one  hand  we  reduce  the  calculation  time  needed  to  solve  the  Volterra 
integral  equation  since  we  are  dealing  with  smaller  distances  in  each  individual  section. 
On  the  other  hand,  it  will  be  necessary  to  calculate  the  field  on  vertical  lines  at 
intermediate  sections.  There  will  be  an  optimum  number  of  sections,  where  the  total 
computation  time  will  be  minimized. 


Figure  2:  Representation  of  the  Multiple  Sections  Solution 
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B.  INITIAL  FIELD  SPECIFICATION 


In  order  to  be  able  to  band  limit  the  initial  data  as  well  as  to  account  for  finite 
antenna  beamwidths  we  use  a  Gaussian  spatial  distribution  given  by: 

(z-zt)2 

u(0,z)  =  Ae  2ff*  (2.7) 


where  z,  is  the  transmitting  antenna  height,  <7Z  is  it  standard  deviation  for  the  Gaussian 
distribution  and  is  related  to  the  antenna  beamwidth  in  the  vertical  plane  through 

BW^  =2  tan-1 

and  A  is  the  amplitude.  Obviously,  choosing  a  small  <7z(for  example  AH)  makes  the 
source  field  very  concentrated  around  z, .  Notice  Figure  37  (upper  figure)  where  we  have 
the  source  field  as  a  function  of  <rz . 

Using  this  source  field  as  f(rj)  in  the  first  integral  in  Equation  2.2  we  find  that  the 
free  space  field,  which  corresponds  to  the  first  integral  of  Equation  2.2,  is  given  by: 


0.833 

h<J, 


M  free- space  Z) 


Act, 


1  IX 


( Z-Ztf 


(2.8) 


Given  the  source  field  described  by  Equation  2.8  we  can  also  calculate  its 
derivative  at  (x,  z)  with  respect  to  z : 


(2.9) 
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C.  NUMERICAL  SOLUTION  PROCEDURE 


1.  Calculation  of  the  Incident  Current  on  the  Rough  Surface 
a.  Calculation  for  the  First  Section 

The  incident  current  on  the  surface  is  given  by  the  first  term  on  the  right 
hand  side  in  Equation  2.6. 


00 

Jinc  fa)  =  -2  J /(?7)—  G0  [x,  g(x);  0,  T]}dT] 


(2.10) 


For  the  first  section  we  have  an  analytical  expression  for  f{r}),  given  by 
Equation  2.7.  Substituting  z  =  g(x) ,  taking  the  derivative  with  respect  to  77  in  Equation 
2.3,  and  using  this  result  and  Equation  2.7  into  the  right  hand  side  of  Equation  2.10,  we 
find  the  expression  for  the  incident  current  on  the  rough  surface.  This  expression  is  valid 
only  for  the  first  section  where  an  analytical  expression  for  the  initial  field  is  available: 


4cfa)= 


-2A<rz[g(*)-zJ 

iK 


2  hr: 


2  « 


2  ** 

°i+T 

Ko, 


(2.11) 


b.  Calculation  for  the  Remainder  Sections 

For  the  remainder  steps  we  need  to  perform  a  numerical  evaluation  of  the 
incident  current,  given  by  Equation  2.10.  The  use  of  trapezoidal  rule  is  not  suitable 
because  the  Green’s  function  term  inside  the  integral  in  Equation  2. 10  has  a  phase  term  of 


the  type  — ,  and  for  small  distances  x  it  introduces  numerical  errors.  As  a  result,  a  more 
x 

sophisticated  integration  technique  is  required. 
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Making  use  of  the  spectral  representation  in  Equation  2.4  we  rewrite  the 
expression  for  Jinc  (x)  as: 

•U*)=— IfWlpe^'e^dp  (2.12) 

k  i 

So 

Interchanging  the  orders  of  integration  and  substituting  p  =  g0+£  we  get: 

ip2x  oo 

ft  0 

=  -t  H  w"f  pe*  3  (p)dp  (2.1 3) 

ft  L 

oo 

where  3(/?)= J  f(g0  +£)e~ip*d£  is  the  Fourier  transform  of  the  initial  field, 
o 

Since  the  parabolic  equation  has  a  limited  angular  response  [Ref.  2],  we 
need  to  truncate  the  integration  limits  in  Equation  2.13  at  ±  ,  where  p^  is  of  the 

order  of  k0sin{6m ax)  where  ^  ^15°-  Because  we  are  truncating  the  integral  in 
Equation  2.13,  we  need  to  introduce  a  filter,  such  as  a  Hanning  filter,  in  order  to  smoothly 
roll-off  the  function.  The  expression  for  J inc  then  becomes: 

Jl„(x)-j-rJw(p)p3(py'^e,'U^ldp  (2.14) 

~P  max 

where  W{p)  is  the  Hanning  filter.  The  above  integral  will  be  evaluated  numerically  by 
the  trapezoidal  rule.  The  Hanning  filter  is  an  even  function  and  its  discrete  form  is  given 
by: 
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W(p  =  mAp)  = 


1,0  <m<-M 
8 

( 4mn')  3 


sin 


N 
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m  <  p  <M 


where  M  is  the  number  of  points  used  in  the  discrete  Fourier  transform.  The  final 
numerical  expression  for  J  inc  is: 


•a  '  M  ! £js± 

i^SW'(p.k3(p.)e 

ft  m-\ 


where  pm=  mAp . 


(2.15) 


2.  Calculation  of  the  Total  Current  on  the  Rough  Surface 

Once  we  calculate  the  incident  current  on  the  rough  surface,  the  next  step  is  to 
calculate  the  total  current  on  the  surface.  The  total  current  on  the  surface  is  represented 

du 

by  Equation  2.6.  We  can  say  for  the  clarity  of  explanation  that  J(x)=-^-(x)  and 

J{£)  =  (£) .  Defining  gd  (*;£)  =  f— ■■  we  can  rewrite  Equation  2.6  as: 

drj  x-g 


J(x)  =  Jinc(x)  + 


'k<)~r~8d(x’4) 


(2.16) 


where  J  inc  (x)  is  given  by  Equation  2.10. 

The  second  term  on  the  right-hand  side  in  Equation  2.16  can  be  broken  in  N  sub¬ 
integrals: 
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def  fit  h 

=  7"c^)+JfV ~lrtn 


(2.17) 


where  /„  =  J 

J  _  J  Y  —  r 


n~  J  yb'  /7Z7 

(n-l)<5  VA  ^ 

Notice  from  Figure  3  below  that  the  following  notation  has  been  adopted: 

n  =  l,2,3v>  A  +  l 

xn  ={n-\)8 


XN+ 1  ~  fit 8  —  x 


x-8  =  xn=(N-1)8 


8n=S(Xn ) 

Jn=J(*n) 


n  1  2 


n  n+1 


x  0  8 


(n-l)5  n8 


Figure  3:  Nomenclature  for  the  discretization  process 
With  reference  to  the  Figure  4,  we  can  approximate  gd{x\%)  = - — — as 

X  £ 


g(x)-g(x,,) 


— ,  where  x  l/=(n-—)8  and  g(xn+w),a  —  Also  we  can 
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approximate  the  variable  £  in  gd(x;£),  located  between  xn  =(n-l)S  and  xn+l  =  nS  by 


(n-l)S  % 


x  =  xN+i 
N8 


Figure  4:  Approximation  for  gd  (x;  t) 

Letting  %  =  T8+{n-\)8 ,  in  Equation  2.17,  where  0<-r<l  in  the  n*  segment, 
and  approximating  /(£)  linearly  between  Jn  and  Jn+l  as  7(|)  =  (/n+1  - )r+ ,  we 
simplify  the  expression  for  I n  as: 


rr  '•y4d,(w+ 1-»)  i 

h=Sdn^ISe1  I 


r-Vr2  r 
1  2 


(2.18) 


The  result  may  be  expressed  in  terms  of  the  Fresnel  integral  as  [Ref.  1]: 
/.  =2^sgn(^){  {((/.., 


(/B+,-7j  a?  r  4*i 
Md  L  J 


for  Sdn* o>  and  /„  =0  for  gdn  =0  (2.19) 


In  Equation  2.22  we  notice  that  I N  =  0  whenever  gN  =  gN+1 . 

Combining  Equations  2.21  and  2.22  we  have  the  final  expression  for  the  current 
on  the  surface: 
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t  _  ^  inc(.N+\) 

J  N+ 1  -  1 -  + 


A 


-2-a,  1 

2m  V  2m 


a, 


l-< 


2m 


.ai 


N-l 


/I=l 


K  Sg 


dn  J 


[Ffc)-Ffe)]+ 


K$sl 


t2e  2  -tx 


(2.23) 


where  the  coefficients  ax  and  bx  are  given  by: 


a,  =2sgn (gN+t-g„)J-^<  F(r2jl+]i 


.JT  2 

^  e'2'2 


^  =2sgn(gw+I-gw). 


i|Z-L. 

V  ^0 


./T  2 


Mil  o 


3.  Calculation  of  the  Field 

Once  we  have  the  total  current  on  the  rough  surface  we  can  calculate  the  field  at 
any  x .  In  Equation  2.2  we  notice  that  w,  (x,  z)  corresponds  to  the  incident  field  and 
u2  (x,  z)  corresponds  to  the  field  reflected  and  diffracted  by  the  rough  surface. 
a.  Calculation  of  the  Incident  Field 

The  incident  field,  corresponding  to  the  first  integral  on  the  right  hand  side 
of  Equation  2.2  and  given  by  ux(x,z)  =  ^f(rj)G0[x,z',O,Tj]dTi  will  be  calculated  using 
the  same  technique  as  used  for  the  calculation  of  J inc  (x).  It  is  given  by: 


ip  x 


ux  (x,  z)  =  —  J 3 (p)e  2k°  eip^z  8o^dp 


(2.24) 


As  before,  a  Hanning  filter  is  introduced  into  the  expression  when  the 
limits  of  the  integral  are  truncated.  The  expression  for  ux{x,z)  then  becomes: 

,  Jil 

ux{x,z)  =  —  f  W{p)Z{p)e  2k°eip[z-g° ]dp  (2.25) 

2  jz  „ 

Pnax 

Finally  we  can  apply  the  trapezoidal  rule  to  Equation  2.25,  obtaining  the  final  numerical 
expression  for  ux(x,z ): 


2^r  „= i 


'p»* 


(2.26) 


b.  Calculation  of  the  Field  Reflected  and  Diffracted  by  the  Rough 


Surface 

The  field  reflected  and  diffracted  by  the  rough  surface,  corresponding  to 
the  second  integral  on  the  right  hand  side  of  Equation  2.2  and  given  by 

Ul(X,Z)=W^l =  (2-27) 

cannot  be  calculated  by  applying  the  trapezoidal  rule  directly  because  the  Green’s 
function  has  a  singularity  at  x  =  % .  This  singularity  must  be  extracted  apriori  before 
performing  the  integration.  We  can  rewrite  the  above  expression  for  u2(x, z)  as: 


r(z\  ,[z-g(£)f  *0 

J_M  -  2(*-#) 


(2.28) 


'8^0JoV^ 

We  replace  [z-g{£))  by  [z-g(x)+g(x:)-g(^)]  in  Equation  2.28  and 


define: 


K{x-,Z)  =  j(Z)e 


»oU-8(*)l 


gfxhsM  ,*ob(*)-g(#)]2 

(*-#)  e  2(x-£) 


(2.29) 
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so  that 


The  integral  in  Equation  2.32  can  now  be  evaluated  by  applying  the 
trapezoidal  rule  to  the  second  term  on  the  right  hand  side,  since  the  singularity  at  x  =  £ 
has  been  extracted. 

4.  Fresnel  Integral  Calculation 

For  small  arguments  (smaller  than  or  equal  than  0.1)  the  Fresnel  integral  can  be 


approximated  by 
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(2.34) 


F(x)=xA  3 
6 

This  expression  is  obtained  by  taking  the  limit  when  t  — >0  in  Equation  2.20.  Using  this 
approximation  will  reduce  CPU  time  somewhat.  For  arguments  greater  than  0.1,  the 
Fresnel  integral  will  be  calculated  by  a  rational  approximation  given  in  [Ref.  3]: 

F(x)  =  C(x)+iS(x) «  —  (l + z')-[z>(*)+,s(x)]e  2  (2.35) 

2 


where 


i*)= 


1+0.926* 


2+1.792*+3.104*J 


■ + f(x) 


?(*)  = 


1 


2+4.142x+3.492*2  +  6.670*; 


■+e(x) 


where  |£(x)|  -  2x10  3 

D.  DETERMINATION  OF  THE  EQUIVALENT  IMPEDANCE 

Once  we  have  determined  the  field  on  a  vertical  line  distant  *  from  the  source, 
we  can  find  the  equivalent  ground  impedance  that  would  produce  the  same  field  when 
solving  the  propagation  problem  using  a  two-ray  model. 

For  horizontal  polarization,  the  reflection  coefficient  for  plane  waves  incident 

77o 

from  air  on  to  a  non-magnetic  lossy,  planar  interface  having  an  impedance  tj=—==  is: 

Verc 

r=sm^-^=j7 ,  (2.36) 

sin^+-^7 


where  if/  is  the  grazing  angle  such  that 


tmy/- 


H,  +Hr 
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:  sin^for  Ht  +Hr  <- 


and  T)0  is  the  intrinsic  impedance  of  the  free-space. 

When  considering  the  two-ray  model  the  resultant  field  will  be  given  by: 

rX  r2  [  h  J 


(2.37) 


where  rx  is  the  line-of-sight  path,  r2  is  the  reflected  path  and  u0  is  the  free-space  field. 

For  large  distances  (d  »  H,  +  Hr )  we  can  approximate  r2  ~  r,  for  the  amplitude 
2  jj 

term,  and  use  r2  -  r{  ~ - 1 — -  for  the  phase  term.  We  can  rewrite  Equation  2.37  as 

d 


u 


un 


=i+rv 


(2.38) 


and  express  the  reflection  coefficient  as: 


r= 


u 

Kuo 


iUpHtHr 


-1 


(2.39) 


If  we  now  combine  Equations  2.36  and  2.39  we  have  the  final  expression  for 
^erc  as  a  function  of  the  geometry  and  the  normalized  field: 


(  nkpH.H,  u  ^ 


^7  =  sin^ 


+  1-- 


un 


1  + 


un 


(2.40) 


Figure  5  below  illustrates  the  geometry  mentioned  above. 
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Figure  5:  Reflection  of  Plane  wave  over  Flat  Surface 
Note  that  the  heights  H,  and  Hr  are  relative  to  the  equivalent  flat  plane,  which  is 

normally  taken  to  coincide  with  the  mean  of  the  rough  surface.  If  only  one  section  is 
used  in  the  PE  modeling  we  can  take  advantage  of  the  Gaussian  representation  of  the 
initial  field.  A  somewhat  better  approximation  to  Equation  2.37  will  then  be: 
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in.  RESULTS 


Having  introduced  the  theoretical  formulation  in  Chapter  IE,  we  will  present  here 
numerical  results  obtained  with  the  algorithm  developed  and  compare  those  with  other 
techniques  for  some  practical  cases  of  radiowave  propagation  over  terrain.  We  will  be 
comparing  the  propagation  factor  (PF)  (excess  loss  over  the  free  space  case)  obtained 
with  our  approach  versus  the  propagation  factor  obtained  from  results  available  in  the 
literature.  A  positive  (negative)  value  of  propagation  factor  implies  gain  (loss)  with 
respect  to  propagation  in  free  space. 

A.  PROPAGATION  OVER  A  PERFECT  ELECTRIC  CONDUCTOR  PLANE 

Consider  a  transmitter  5  m  above  flat  perfect  electric  conductor  (PEC)  ground 
(z,  =5).  One  wants  to  calculate  the  PF  on  a  vertical  line  distant  400  m  from  the 
transmitter,  as  shown  in  Figure  6.  The  maximum  receiver  height  to  be  evaluated  is  100 
m. 

!  i  7T 

100  m 


(3.1) 


zt  =  5  m 


x  =  400m 


dl 


-Xr 


d2 


dn 

< - > 


< - 

Figure  6:  Propagation  over  a  Flat  PEC  Surface 
The  PF  is  defined  as: 

f  m(x,z)  ^ 


PF  [<®]=101ogIO 


u0{x,z) 


H  = 
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where  u{x,z)  and  u0  (x,  z)  are  defined  in  Equations  2.2  and  2.8  respectively. 

For  a  point  source  above  a  PEC  plane  an  analytical  solution  can  be  obtained  by 
image  theory.  Using  the  expression  valid  in  the  far  field  for  x  »  z, ,  Z ,  we  get 


u 

f  hzz,' 

=  2 

sin 

u0 

{  x  \ 

Figure  7  shows  a  comparison  between  the  two  methods  for  a  frequency  of 
operation  of  1  GHz  and  performing  the  computation  using  one  section.  For  this  range  of 
400  m  good  agreement  is  seen  for  heights  up  to  60  m.  This  result  is  expected,  since  the 
2-ray  model  and  the  PE  model  use  approximations  that  are  valid  for  horizontal  distance 
large  compared  to  antenna  heights.  If  we  move  the  receiver  vertical  line  further  to  the 
right  on  Figure  6,  thereby,  increasing  the  distance  from  the  source,  we  will  have  better 
agreement  at  higher  heights. 

On  Table  3.1  we  have  the  relationship  between  the  CPU  time  (Pentium  MMX  200 
MHz)  and  the  number  of  sections  used  in  the  calculation  of  the  field  over  a  total  distance 
of  400  m.  The  code  was  implemented  in  MATLAB.  The  column  %  time  refers  to  a 
relative  measure  to  the  maximum  measured  time.  We  notice  that  as  we  change  the 
number  of  sections  from  one  to  two  and  from  two  to  three  the  CPU  time  decreases. 
However,  as  we  move  from  three  to  four  sections  the  CPU  time  increases.  Of  course 
changing  the  section  length  will  make  the  CPU  time  change  too.  Compare  in  Table  3.1 
cases  two  and  three,  both  using  2  sections;  cases  four  and  five,  both  using  3  sections  and 
cases  six  and  seven,  both  using  4  sections.  Figure  8  illustrates  this  dependence.  These 
numbers  are  influenced  by  many  factors,  some  dependent  on  the  set-up  of  the  problem 
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(number  of  points  on  the  vertical  direction,  for  example)  and  some  dependent  on  the 
computer  configuration  (available  memory,  for  example). 

We  also  notice  that  it  is  not  possible  to  choose  the  sections  too  small,  since  it 
would  mean  that  waves  would  have  to  travel  at  large  angles  to  reach  large  vertical  heights 
and  the  parabolic  equation  has  an  angular  limitation  around  15°.  Referring  to  Figure  9 
we  notice  that  the  results  for  the  cases  3,  5  and  7  are  smoother  than  2,  4  and  6 
respectively.  This  is  because  the  fields  on  a  vertical  line  get  more  accurate  as  the 
horizontal  distance  from  the  source  increases. 

We  can  make  a  better  check  to  our  numerical  solution  by  comparing  it  with  the 
solution  for  Gaussian  sources  instead  of  the  point  sources  used  above.  For  a  Gaussian 
source,  described  by  Equation  2.7,  an  analytical  expression  for  the  total  field  over  a 
ground  plane  is: 


where  the  first  term  on  the  right  hand  side  corresponds  to  the  free-space  field  and  the 
second  term  corresponds  to  the  field  generated  by  the  image  source.  The  PF  will  be  given 


Figure  10  shows  current  and  the  propagation  factor  for  Case  2  (dl  =  d2  =  200  m). 
There  we  notice  a  small  oscillatory  behavior  for  the  current.  Compare  this  figure  with 
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Figure  7,  where  the  calculations  were  performed  in  one  single  section.  It  is  seen  that  the 
results  are  smoother.  Below  we  offer  an  analysis  of  this  behavior. 

Figure  1 1  compares  the  propagation  factor  at  the  end  of  the  first  section  in  Case  2 
(dl  =  d2  =  200  m),  when  using  our  numerical  solution  and  using  the  image  theory  with 
Gaussian  source  as  described  by  Equation  3.2.  We  notice  such  a  good  agreement,  even 
for  large  heights,  that  it  is  hard  to  distinguish  the  two  curves.  Figure  12  compares  the 
magnitude  and  the  phase  of  the  field  at  the  end  of  the  first  section.  For  the  sake  of  clarity 
we  present  Figure  13  that  contains  basically  the  same  information  as  Figure  12  but  over  a 
smaller  vertical  scale.  There  we  compare  our  numerical  solution  with  the  solution 
obtained  using  the  image  theory  described  above.  We  notice  a  very  good  agreement, 
which  indicates  that  the  field  at  the  end  of  the  first  section  is  being  calculated  correctly. 
As  a  result,  we  conclude  that  although  the  multiple  section  solution  improves  the  speed  of 
the  calculation,  it  introduces  some  numerical  errors,  since  the  PF  for  Case  1  and  Case  2 
are  not  the  same,  and  the  only  difference  between  the  two  cases  is  that  in  Case  2  we 
performed  the  computation  in  two  sections.  In  Section  3.C  we  will  see  that  as  we 
decrease  the  length  of  the  first  section  the  results  become  worse. 
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Propagation  Factor  (dB) 


Figure  7:  Current  (Upper  Figure)  and  Propagation  Factor  (Lower  Figure).  Transmitter 

A 

antenna  height  is  5  m.  Flat  earth,  frequency  =  1  GHz,  A  =  — . 


CPU  Time  (Seconds) 


Case 


N.  of 

sections 


Section  Size  (m) 


CPU  time 


1.0177  x  10J 
743.5800 
833.7700 
638.2300 
544.2000 
646.9100 
651.4100 


Table  3.1:  CPU  time  for  Propagation  over  PEC  Flat  Surface 


Case  1: 
Case  2: 
Case  3: 
Case  4: 
Case  5: 
Case  6: 
Case  7: 


1  section 

2  sections 

2  sections 

3  sections 

3  sections 

4  sections 
4  sections 


2  2.5  3 

Number  of  Sections 


Figure  8:  CPU  Time  as  a  Function  of  Number  of  Sections  for  the  PEC  Flat  Surface  Case 


^  |  <N 


Receiver  Height  (m) 


Figure  11:  Propagation  Factor,  Calculated  at  the  End  of  the  First  Section  of  Case  2  of  Flat 


Surface. 


Figure  13:  Magnitude  and  Phase  of  the  Field  at  the  End  of  the  First  Section  for  Case  2  of 
Hat  Surface. 
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B.  PROPAGATION  OVER  PEC  KNIFE-EDGE 


Consider  an  ideal  absorbing  knife-edge  of  height  hp  =  5  meters  on  a  conducting 

plane,  as  shown  in  Figure  14,  where  we  have  a  =  b  =  c  =  d  =  50m.  The  transmitter 
antenna  height  is  3  meters  ( z,  =  3 )  and  it  is  located  at  point  A.  One  wants  to  calculate 
the  PF  at  a  vertical  line  located  at  B,  distant  200  m  from  the  transmitter  antenna. 


Figure  14:  Perfectly  Absorber  Knife-Edge  Between  the  Transmitter  at  A  and  a  Vertical 
Line  at  B,  Where  the  Propagation  Factor  is  Being  Evaluated.  Both  are  over  Perfectly 
Conducting  Ground. 

We  compare  our  numerical  solution  with  a  four-ray  model  of  knife-edge 

diffraction  theory  [Ref.  4].  In  the  parabolic  equation  computations  we  replace  the  knife- 

edge  with  a  triangular  hill  having  a  base  of  100  m  (b  =  c  =  50  m  in  Figure  14)  and  a  peak 

height  of  5  m.  The  numerical  results  are  plotted  in  Figure  15.  We  notice  that  the  total 

current  and  the  incident  current  are  the  same  for  the  first  50  m  (corresponding  to  portion 

a),  since  we  have  flat  surface  for  this  part  of  the  terrain.  We  notice  a  good  agreement  for 

low  heights,  since  the  four-ray  model  assumes  large  distances.  If  we  now  move  both  the 

triangular  peak  and  the  vertical  line  B  to  the  right,  making  a  =  d  =  150  m  and  b  =  c  =  50 
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m  on  Figure  14  we  have  the  numerical  results  shown  in  Figure  16.  Moving  the  triangular 
peak  and  the  vertical  line  B  further  to  the  right,  making  a  =  d  =  250  m  and  b  =  c  =  50  m 
on  Figure  14,  we  have  the  results  shown  on  Figure  17.  Comparing  Figure  15,  Figure  16 
and  Figure  17  we  notice  progressively  better  agreement  with  increasing  heights,  since 
both  models  compared  provide  more  accurate  results  for  longer  horizontal  distances. 
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Figure  15:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  figure)  for  the  Triangular  Hill  with  a  =  b  =  c  =  d  =  50m,  Peak  Height  =  5  m  and 
Antenna  Height  =  3  m.  Frequency  =  1  GHz.  Calculation  in  one  section.  Ax  =  A . 
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2.5 


-  Total  Current 


Figure  16:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  figure)  for  the  Triangular  Hill  with  a  =  d  =  150  m,  b  =  c  =  50  m,  Peak  Height  =  5 
m  and  Antenna  Height  =  3  m.  Frequency  =  1  GHz.  Calculation  in  one  section.  Ax=X. 


35 


Figure  17:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  figure)  for  the  Triangular  Hill  with  a  =  d  =  250  m,  b  =  c  =  50  m,  Peak  Height  =  5 
m  and  Antenna  Height  =  3  m.  Frequency  =  1  GHz.  Calculation  in  one  section.  Ax=X. 
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C.  PROPAGATION  OVER  A  PEC  SINUSOIDAL  SURFACE  AND  EQUIVALENT 
IMPEDANCE  DETERMINATION 

Once  we  tested  our  solution  for  two  important  test  cases,  we  can  use  the  model  for 
a  rough  surface.  The  only  restriction  we  enforce  is  on  the  slope  of  the  surface:  as  it 
increases,  the  waves  are  forced  to  travel  in  higher  angles  and  the  parabolic  equation  has 
an  angular  limitation  of  approximated  15°  [Ref.  2].  In  this  section  we  will  analyze  the 
case  of  a  sinusoidal  surface.  Initially  we  will  evaluate  the  propagation  factor  and  verify 
how  fast  we  can  march,  i.e.,  how  the  results  are  sensitive  to  changes  of  the  step  size. 
Finally,  we  will  determine  the  equivalent  impedance  for  low  grazing  angles. 

1.  Propagation  over  a  PEC  sinusoidal  surface 
Consider  a  sinusoidal  surface  defined  by 

(  \  App(,  f  2;zx 

g{x)  =  ^T-  l"cos  — +e 
1  V  v 1 

where  App  is  the  peak  to  peak  amplitude,  x  is  the  horizontal  distance  from  the  source, 
T  is  the  period,  6  an  arbitrary  phase  and  g(x)  the  amplitude  on  that  particular  horizontal 
distance. 

In  Table  3.2  we  have  the  relationship  between  the  CPU  time  (Pentium  MMX  266 
MHz,  MATLAB  code)  and  the  number  of  sections  used  in  the  calculation  of  the  field  over 
a  total  distance  of  800  m  on  the  sinusoidal  surface  with  the  following  parameters: 
A  =  lm,  T  =  30  m,  0  =  0.  The  distances  dl,  d2,  d3  and  so  on  are  defined  in  Figure  18. 
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Table  3.2:  CPU  time  for  Propagation  over  PEC  Sinusoidal  Surface 


Figure  19  shows  the  incident  current,  the  total  current  and  the  propagation  factor 
when  performing  the  calculation  in  one  section.  Figure  20,  Figure  21  and  Figure  22  show 
the  incident  current,  the  total  current  and  the  propagation  factor  when  performing  the 
calculation  in  two  sections,  gradually  decreasing  the  length  of  the  first  section  (dl)  and 
increasing  the  length  of  the  second  section  (d2).  We  notice  that  as  we  decrease  the  length 
of  the  first  section  and  increase  the  length  of  the  second  section  we  have  a  disagreement 
among  the  results.  Comparing  cases  2  and  3  with  case  1  the  disagreement  is  not 
significant.  However,  case  4  and  case  1  have  completely  different  results. 
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Figure  23  shows  the  incident  current,  the  total  current  and  the  propagation  factor 
in  three  sections.  Figure  24  shows  the  incident  current,  the  total  current  and  the 
propagation  factor  in  nine  sections. 

Comparing  Figure  19  and  Figure  24,  where  we  have  the  calculation  in  one  section 
and  nine  sections  respectively,  we  notice  that  the  currents  and  the  PF  become  very  noisy 
for  the  last  case.  However,  we  have  the  computation  time  reduced  to  41  %  of  the  original 
one.  In  contrast  to  the  case  of  flat  surface,  we  observe  that  for  propagation  over  rough 
surface  we  need  to  have  many  small  sections  to  reach  the  point  such  that  the 
computational  time  becomes  minimum.  This  does  not  happen  for  the  flat  surface  case 
because  we  had  all  g  dr  ’s  equal  to  zero  and  we  did  not  have  to  calculate  the  Fresnel 

integral  when  evaluating  the  I n  ’s  (refer  to  Equation  2.19). 
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Figure  19:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  Figure)  for  the  Case  1  of  Sinusoidal  Surface.  Antenna  Height  =  5  m.  Frequency 
=  1  GHz. 
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Distance  (m) 


Figure  20:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  Figure)  for  the  Case  2  of  Sinusoidal  Surface.  Antenna  Height  =  5  m.  Frequency 
[  GHz. 
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(Lower  Figure)  for  the  Case  3  of  Sinusoidal  Surface.  Antenna  Height  =  5  m.  Frequency 
=  1  GHz. 
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Propagation  Factor  (dB) 


Figure  22:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  Figure)  for  the  Case  4  of  Sinusoidal  Surface.  Antenna  Height  =  5  m.  Frequency 


Propagation  Factor  (dB) 


Figure  23:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  Figure)  for  the  Case  5  of  Sinusoidal  Surface.  Antenna  Height  =  5  m.  Frequency 


Figure  24:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  Figure)  for  the  Case  10  of  Sinusoidal  Surface.  Antenna  Height  =  5  m.  Frequency 
=  1  GHz 
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2.  Sensitivity  to  the  step  size,  A,,  and  to  the  standard  deviation, cz,  of  the 
Gaussian  source  Held 

Consider  propagation  over  a  sinusoidal  surface  with  the  following  parameters: 
zt  =  5  m ,  App  =  1  m ,  T  =  30  m ,  x  =  300 m ,  6  =  0°  and  frequency  of  1  GHz  ( A  =  0.3 m). 

Figure  25,  Figure  26,  Figure  27,  Figure  28,  Figure  29  and  Figure  30  shows  the  incident 

A 

current,  the  total  current  and  the  PF  for  Ax  equal  to  —  (0.15  m),  A  (0.3  m),  2 A  (0.6  m), 

2 


4 A  (l.2m),  8 A  (2.4  m)  and  1 6 A  (4.8  m)  respectively  and  the  source  field  with  crz  =0.1. 
Figure  31,  Figure  32,  Figure  33,  Figure  34,  Figure  35  and  Figure  36  shows  the  incident 

A 

current,  the  total  current  and  the  PF  for  Ax  equal  to  —  (0.15 m),  A  (0.3m),  2 A  (0.6m), 
4  A  (l.2m),  8 A  (2.4  m)  and  16 A  (4.8  m)  respectively  and  the  source  field  with 


<7Z  =0.4. 


We  notice  that  for  crz  =0.1  we  have  good  results  for  both  the  total  current  and 
the  PF  for  Ax  up  to  2 A .  For  az  =  0.4  we  have  good  results  for  the  total  current  for 
A ^ up  to  16/1 ,  although  even  for  small  Ax  {A/ 2)  we  have  problems  with  the  PF  for  high 


heights.  This  happens  because  the  free-space  field  becomes  very  small  for  those  heights 
due  to  the  Gaussian  distribution  and  numerical  errors  are  introduced.  In  Figure  37  we 
notice  that  for  oz  =0.1m  the  free  space  field  has  almost  constant  magnitude  at  the 
receiver  location,  while  for  <jz  =  0.4  m  the  field  magnitude  decays  as  we  go  higher.  The 
field  becomes  too  small  for  high  heights  and  numerical  errors  are  introduced  when 
computing  the  PF. 


47 


Magnitude 


Figure  25:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

X  X 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax=~,  <JZ=— 

Jm4  J 
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Figure  26:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

A 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  =  A ,  <TZ  =  — 
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Propagation  Factor  (dB) 


Figure  27:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

X 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax  =21,  crz  =  — 


Fbceiver  Height  (rri)  IVtegptude 
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Figure  28:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  AX=4Z,  <TZ  = 
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|  <r> 


Receiver  Heic^t  (m)  IVbgiitude 


Propagation  Factor  (dB) 

Figure  29:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

X 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax  =  8/1 ,  o z  =  — 


Receiver  Height  (rri)  IVbgitude 


Figure  30:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

A 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax  =  16 A ,  oz  =  — 
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Propagation  Factor  (dB) 


Figure  31:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 
(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax  =— ,  crz  =  — 


Receiver  Height  (nil  IVbyitude 


Figure  32:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

,  4A 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  A x=  A,  (Tz-  — 
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Figure  33:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

4A 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax  =  2X ,  c z  — 
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Propagation  Factor  (dB) 


Figure  34:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

4/1 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax  =  4/L ,  =  — 


Receiver  Hei^t  (m)  Magitude 


Propagation  Factor  (dB) 


Figure  35:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 

42 

(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  Ax  =82,  <7.  =  — 
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Fteceiver  Heicft  (m)  IVb^itude 


Figure  36:  Total  Current  and  Incident  Current  (Upper  Figure)  and  Propagation  Factor 


(Lower  Figure)  for  Propagation  over  Sinusoidal  Surface.  A  x  =  16/1 ,  <TZ  =  — 


Fleeter  He®*  (rr&  Source  Hei^t  (m) 


Free  Space  Field 

Figure  37:  Normalized  Source  Field  (Upper  Figure)  and  Free  Space  Field  (Lower  Figure) 
at  x  =  300  m .  Transmitter  Antenna  Height  =  5  m.  Frequency  =  1  GHz 


3.  Equivalent  Impedance  Determination 

Once  we  have  calculated  the  propagation  factor  for  the  sinusoidal  surface,  we  can 
find  the  equivalent  impedance  of  the  sinusoidal  surface  for  low  grazing  angles.  If  we 
compare  Figure  7  and  Figure  19  we  notice  that  for  small  receiver  heights  (=  low  grazing 
angles)  both  curves  have  approximately  the  same  shape.  This  can  give  us  an  intuitive 
explanation  that  under  these  circumstances  we  can  represent  the  PEC  sinusoidal  surface 
by  the  lossy  flat  surface. 

Figure  38  shows  the  real  part  and  the  imaginary  part  of  the  normalized  equivalent 


impedance 


for  the  sinusoidal  surface  using  the  following  parameters:  z,  =  1.5  m , 


A  =  1  m ,  T  =  30  m  and  x  =  300  m .  We  ran  the  problem  50  times  with  6  uniformly 

distributed  in  the  interval  [0,2;r] .  We  can  see  that  as  the  grazing  angle  approaches  zero, 
the  real  part  tends  to  1.69  and  the  imaginary  part  tends  to  4.58.  Notice  the  dashed  line  in 
the  first  plot  on  Figure  38,  which  represents  the  best  linear  fit  when  minimizing  the  sum 
of  the  squared  error.  If  we  now  run  the  flat  surface  problem  with  these  values  we  obtain 
the  same  results  for  PF  as  the  sinusoidal  surface  when  averaged  over  many  realizations. 
Notice  the  bottom  plot  on  Figure  38.  Figure  39  shows  results  when  running  100 
realizations.  We  notice  that  as  the  grazing  angle  approaches  zero,  the  real  part  tends  to 
2.16  and  the  imaginary  part  tends  to  4.56. 

In  Figure  40  we  ran  50  realizations  with  the  parameters  described  above 
(z,  =  1.5  m ,  App=  1  m ,  T  -  30  m ,  x  =  300  m  and  6  uniformly  distributed  in  the  interval 

[0,2#]),  but  we  changed  the  maximum  receiver  height  to  78  m,  or,  equivalently,  the 

maximum  grazing  angle  to  15°.  Then,  using  the  equivalent  impedance  found  when 
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performing  the  fitting  for  low  grazing  angles  (real  part  tending  to  1.69  and  imaginary  part 
tending  to  4.58)  we  calculated  the  equivalent  PF  for  propagation  over  flat  lossy  surface. 
In  the  upper  figure  we  notice  a  linear  behavior  for  both  the  real  and  imaginary  part  of  the 
equivalent  impedance  for  a  maximum  grazing  angle  of  approximately  7°.  In  the  lower 
figure  we  notice  a  good  agreement  between  the  propagation  factor  when  performing  50 
realizations  and  the  equivalent  flat  lossy  surface  for  a  maximum  height  of  approximately 
15  m.  This  result  is  expected,  since  as  we  observed  in  Figure  7  and  Figure  19,  we  have 
approximately  the  same  shape  only  for  small  receiver  heights.  In  other  words,  the 
substitution  of  the  sinusoidal  PEC  surface  for  a  lossy  flat  surface,  when  comparing  the 
PF,  is  valid  only  for  low  grazing  angles. 

For  extremely  small  grazing  angles,  the  errors  in  —  are  amplified  by  Equation 

Uq 

2  39,  due  to  the  presence  of  — — .  This  is  the  cause  of  the  large  swings  in  Figure  38, 

sin^ 

Figure  39,  and  Figure  40  (upper  figures)  at  low  grazing  angles  and  is  also  the  cause  of 
discrepancy  seen  in  the  PF  in  the  PF  in  Figure  38,  Figure  39,  and  Figure  40  (lower 
figures).  Notice  that  we  have  smaller  swings  for  one  hundred  realizations  than  fifty 
realizations. 
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Figure  38:  Upper  Figure:  Equivalent  Impedance  of  Sinusoidal  Surface  with  Period  =  30 
m  and  Peak  to  Peak  Amplitude  =  1  m.  Total  Distance  =  300  m.  Frequency  =  1  GHz. 
Antenna  Height  =  1 .5  m.  Ax  =  0.3m.  Fifty  realizations.  Lower  Figure:  PF  for  the  PEC 
Sinusoidal  Surface  Above  and  for  the  Equivalent  Flat  Lossy  Surface. 
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Figure  39:  Upper  Figure:  Equivalent  Impedance  of  Sinusoidal  Surface  with  Period  =  30 
m  and  Peak  to  Peak  Amplitude  =  1  m.  Total  Distance  =  300  m.  Frequency  =  1  GHz. 
Antenna  Height  =  1.5  m.  Ax  =  0.3m.  One  hundred  realizations.  Lower  Figure:  PF  for 
the  PEC  Sinusoidal  Surface  Above  and  for  the  Equivalent  Flat  Lossy  Surface. 
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Figure  40:  Upper  Figure:  Equivalent  Impedance  of  Sinusoidal  Surface  with  Period  =  30 
m  and  Peak  to  Peak  Amplitude  =  1  m.  Total  Distance  =  300  m.  Frequency  =  1  GHz. 
Antenna  Height  =  1.5  m.  Ax  =  0.3m.  Fifty  realizations.  Lower  Figure:  PF  for  the  PEC 
Sinusoidal  Surface  Above  and  for  the  Equivalent  Flat  Lossy  Surface. 
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IV.  SEA  SURFACE  GENERATION 


The  sea  spectrum  can  be  modeled  in  several  different  ways  depending  on  the 
parameters  that  one  tries  to  represent.  In  this  chapter  we  will  use  the  Pierson-Moskowitz 
spectrum,  which  represents  the  sea  spectrum  as  a  function  of  the  wind  speed  and  is  given 
by  [Ref.  5]: 


k2U< 


(4.1) 


where  a  =  8.1xlO~3 ,  ft  =  0.74,  g  =  9.81  m/s2 ,  U  =wind  speed  (m/s)  at  a  height  of  19.5 
2k 

m,  k  —  —  is  the  wave  number  and  A  is  the  surface  wavelength. 

A 

Figure  41  shows  the  spectrum  for  wind  speed  of  5  m/s  and  6  m/s . 


01  23456789 

Wave  Number  (1/m) 


Figure  41:  Pierson-Moskowitz  Spectrum  for  Wind  Speed  of  5 m/s  and  6m/s 
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Realizations  of  the  random  surface  will  be  generated  at  points 
xn  =  nAx  (n  =  1,...,  N)  according  to  [Ref.  6]: 

/k)=7  <4-2) 

where,  for  j  >  0 , 

Tw(o.i)+w((u)]_  » 

2*£W(j:)  ^  2  (4.3) 

MOA  y=o,| 

and,  for  7  <  0 ,  f(k,)=  F’(Kj).  In  Equation  4.3  Kt=^y-,  where  L  is  the  surface 

length  and  each  time  iV(0,l)  appears  it  represents  an  independent  sample  taken  from  a 
zero  mean,  unit  variance  Gaussian  distribution.  Equation  4.2  is  computed  with  a  fast 
Fourier  transform. 

Figure  42  (upper  figure)  shows  a  possible  surface  realization  for  a  wind  speed  of 
5  m/s .  Since  the  parabolic  equation  has  a  limited  angular  response  [Ref.  2]  we  can  avoid 
high  slopes  smoothing  the  surface,  i.e.,  changing  a  local  point  by  the  average  between 
itself  and  some  neighbors  points.  Figure  42  (lower  figure)  shows  the  realization  of  the 
upper  figure  when  replacing  any  point  by  the  moving  average  of  itself  and  four 
neighbors,  two  to  the  right  and  two  to  the  left.  When  making  use  of  this  smoothing 
process  we  practically  do  not  change  the  statistics  of  the  surface.  Figure  43  (upper 
figure)  shows  a  comparison  between  the  inverse  Fourier  transform  of  the  noisy  spectrum 
and  the  autocorrelation  of  the  rough  surface  for  the  wind  speed  of  5  m/s .  Figure  43 
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(lower  figure)  shows  the  same  comparison,  but  replacing  the  autocorrelation  of  the  rough 
surface  by  the  autocorrelation  of  the  smoothed  rough  surface. 

Defining  the  correlation  length  7’  as  the  length  such  that  0.5  xZ  equals  the  area 
under  the  autocorrelation  function  (only  for  positive  lags)  and  dx  as  the  separation 
between  two  consecutive  points  on  the  rough  surface,  we  need  to  have 

L»l  and  l»dx 

in  order  to  have  a  truly  uncorrelated  surface.  In  the  example  above,  where  we  had  the 
wind  speed  of  5 m/s ,  we  had  L  -  300  m ,  /  =  21 .2  m  and  dx  =  0.29  m . 

Figure  44  shows  a  surface  realization  (top  figure)  and  the  respectively  smoothed 
surface  realization  (bottom  figure)  for  a  wind  speed  of  20 m/s.  Notice  that  the  surface 
heights  increase  substantially  when  compared  to  Figure  42.  Figure  45  (upper  figure) 
shows  a  comparison  between  the  inverse  Fourier  transform  of  the  noisy  spectrum  and  the 
autocorrelation  of  the  rough  surface  for  the  wind  speed  of  20 m/s.  Figure  45  (lower 
figure)  shows  the  same  comparison,  but  replacing  the  autocorrelation  of  the  rough  surface 
by  the  autocorrelation  of  the  smoothed  rough  surface.  For  this  case  we  of  wind  speed  of 
20m/s  we  have  L  =  2000  m ,  l  =  469.9  m  and  dx  =  3.9  m . 
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Horizontal  Distance  (m) 

Figure  42:  Sea  Surface  Realization  for  Wind  Speed  of  5m/s  (Upper  Figure)  and 
Smoothed  Sea  Surface  (Averaging  Five  Points)  for  Wind  Speed  of  5m/s  (Lower  Figure). 
Length  =  300m 
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Figure  43:  Normalized  Fourier  Transform  of  the  Noisy  Spectrum  and  Rough  Surface 
Autocorrelation  (Upper  Figure)  and  Normalized  Fourier  Transform  of  the  Noisy 
Spectrum  and  Smoothed  Rough  Surface  Autocorrelation  (Lower  Figure).  Wind  Speed  = 
5m/s 
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Figure  44:  Sea  Surface  Realization  for  Wind  Speed  of  20  m/s  (Upper  Figure)  and 
Smoothed  Sea  Surface  (Averaging  Five  Points)  for  Wind  Speed  of  20  m/s  (Lower 
Figure).  Length  =  4000  m 
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Figure  45:  Normalized  Fourier  Transform  of  the  Noisy  Spectrum  and  Rough  Surface 
Autocorrelation  (Upper  Figure)  and  Normalized  Fourier  Transform  of  the  Noisy 
Spectrum  and  Smoothed  Rough  Surface  Autocorrelation  (Lower  Figure).  Wind  Speed  = 
20  m/s 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  thesis,  a  numerical  method  to  model  propagation  over  random  rough 
surface  using  the  parabolic  equation  was  implemented  and  tested.  The  parabolic 
equation,  where  the  Helmholtz  equation  is  approximated  by  a  parabolic  equation,  is  a 
quasi-full  wave  method  and  takes  into  account  aspects  such  as  forward  reflection, 
diffraction,  refraction  and  surface  wave  propagation.  However,  it  ignores  back- 
scattering.  Because  of  this  assumption,  it  is  possible  to  use  a  marching  scheme,  where 
the  field  at  a  certain  point  is  calculated  uniquely  based  on  the  field  at  a  previous  point, 
resulting  in  a  very  efficient  computational  method.  Besides,  this  method  seems  to  be  the 
only  practical  method  for  predicting  radiowave  propagation  over  long  ranges  (thousands 
of  wavelengths)  and  over  a  wideband  of  frequencies. 

Since  our  technique  uses  a  Volterra  integral  equation  of  the  second  kind,  the 
numerical  computations  become  very  intense  for  large  distances.  In  order  to  minimize 
the  computation  time,  we  introduced  a  multiple  sections  solution,  where  the  numerical 
solution  scheme  is  successively  repeated  for  smaller  ranges.  Given  the  field  on  a  vertical 
line  between  two  sections  we  calculate  the  current  on  the  surface  between  these  two 
vertical  lines  and  then  calculate  the  field  on  the  next  vertical  line.  When  using  this 
multiple  sections  solution  we  substantially  decrease  the  computation  time,  although  we 
noticed  that  when  sections  become  too  small  numerical  errors  were  introduced. 

We  tested  our  numerical  techniques  for  some  test  cases  encountered  in  the 
literature  such  as  propagation  over  a  PEC  flat  surface  and  over  an  ideal  absorbing  knife- 
edge.  The  agreement  between  the  numerical  solution  and  these  test  cases  was  good. 


75 


Next  we  tested  our  numerical  solution  for  propagation  over  a  sinusoidal  surface. 
We  noticed  that  for  low  grazing  angles  the  propagation  factor  has  a  shape  similar  to  the 
one  when  propagation  takes  place  over  a  flat  lossy  surface.  This  fact  served  as 
motivation  for  modeling  the  difficult  problem  of  propagation  over  a  rough  surface  at  low 
grazing  angles  by  propagation  over  an  equivalent  impedance  plane.  We  simulated  a 
random  surface  with  a  fixed  spectrum  by  translating  the  sinusoidal  surface  of  given  peak 
to  peak  height  and  period  in  a  random  fashion,  with  a  uniform  distribution  in  [0,2k]. 

Then  we  compared  the  resulting  normalized  field  to  the  field  of  a  2-ray  model  over  a 
constant  impedance  plane.  For  low  grazing  angles  the  agreement  was  excellent.  Finally 
we  showed  how  to  generate  rough  sea  surfaces  according  to  the  Pierson-Moskowitz 
spectrum,  which  represents  the  sea  spectrum  as  a  function  of  the  wind  speed. 

This  thesis  represents  an  effort  to  predict  radiowave  propagation  over  a  random 
rough  surface  such  as  the  sea.  The  next  step  should  be  to  predict  the  propagation  factor 
for  propagation  over  the  sea  surface.  Finally,  because  the  solution  of  the  Volterra  integral 
equation  takes  significant  amount  of  CPU  time  during  recursions,  the  source  code  should 
be  implemented  in  a  programming  language  that  deals  with  for  loops  efficiently,  such  as 
C  or  FORTRAN.  In  this  thesis  we  used  MATLAB,  which  is  known  to  be  very  efficient  for 
matrix  manipulation,  but  very  inefficient  for  loops. 
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